
##############################################
########### IV MODELS w/ DK SE's  ############
##############################################
## Code below produces Appendix Figures AF3 and AF4, and Tables AT2 and AT3

bcs.dat <- read.csv("....R_models_IV.csv")
colnames(bcs.dat)

library(plm)
library(stargazer)
library(ggplot2)
library(egg)
library(ggpubr)

#### Marginal effects plot function, to incorporate the Driscoll-Kraay standard errors ####
dhplot <- function(model,var1,var2,int,vcov,ci=.95,
                   xlab=var2,ylab=paste("Marginal Effect of",var1),
                   main="",
                   me_lty=1,me_lwd=1,me_col="navyblue",
                   ci_lty=1,ci_lwd=.5,ci_col="navyblue",
                   yint_lty=1,yint_lwd=1,yint_col="red"){
  require(ggplot2)
  alpha <- 1-ci
  z <- qnorm(1-alpha/2)
  beta.hat <- coef(model)
  cov <- vcov
  z0 <- seq(min(model.frame(model)[,var2],na.rm=T),max(model.frame(model)[,var2],na.rm=T),length.out=1000)
  dy.dx <- beta.hat[var1] + beta.hat[int]*z0
  se.dy.dx <- sqrt(cov[var1,var1] + z0^2*cov[nrow(cov),ncol(cov)] + 2*z0*cov[var1,ncol(cov)])
  upr <- dy.dx + z*se.dy.dx
  lwr <- dy.dx - z*se.dy.dx
  ggplot(data=NULL,aes(x=z0, y=dy.dx)) +
    labs(x=xlab,y=ylab,title=main) +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_line(aes(z0, dy.dx),size = me_lwd, 
              linetype = me_lty, 
              color = me_col) + 
    geom_ribbon(aes(ymin=upr,ymax=lwr), 
                fill="skyblue4", alpha=0.5) +
    geom_line(aes(z0, lwr), size = ci_lwd, 
              linetype = ci_lty, 
              color = ci_col) +
    geom_line(aes(z0, upr), size = ci_lwd, 
              linetype = ci_lty, 
              color = ci_col) +
    geom_hline(yintercept=0,linetype=yint_lty,
               size=yint_lwd,
               color=yint_col, alpha=0.5)
}

##########################################
############### CBI IV MODELS ############
##########################################

bc.full <- pdata.frame(bcs.dat, index = c("country", "year"))

####
#### UNEMP 
####

iv.cbi1 <- plm(D_logunemp_ilo ~ LD_logunemp_ilo + BCsys * cbi + 
                 D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
               + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis  | LD_logunemp_ilo +
                 D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
               + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis   
               +  L_BCsys * L_CBI_IV_PCA  , 
               data = bc.full, model = "random" )

covDK.un <- vcovSCC(iv.cbi1)
covDK2.un <- sqrt(diag(covDK.un))
stargazer(iv.cbi1,  se=list(covDK2.un), type = "text", star.cutoffs = c(0.1, 0.05, 0.001), star.char = c("*", "**", "***")) 
cbi.Un <-  dhplot(model = iv.cbi1, var1 = "BCsys", var2 = "cbi", int = "BCsys:cbi", vcov = covDK.un, ci=0.90,
                  main = "(D) Unemployment" , ylab = "(Change) in Unemployment", xlab = "CBI",yint_lwd = 0.7)

####
#### CREDIT 
####

iv.cbi2 <- plm(D_logdomcredit ~ LD_logdomcredit + BCsys * cbi + 
                 D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
               + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis  | LD_logdomcredit +
                 D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
               + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis   
               +  L_BCsys * L_CBI_IV_PCA  , 
               data = bc.full, model = "random" )


covDK.cr <- vcovSCC(iv.cbi2)
covDK2.cr <- sqrt(diag(covDK.cr))
stargazer(iv.cbi2,  se=list(covDK2.cr), type = "text", star.cutoffs = c(0.1, 0.05, 0.001), star.char = c("*", "**", "***")) 
cbi.cr <-  dhplot(model = iv.cbi2, var1 = "BCsys", var2 = "cbi", int = "BCsys:cbi", vcov = covDK.cr, ci=0.90,
                  main = "(D) Credit/GDP" , ylab = "(Change) in Credit/GDP", xlab = "CBI",yint_lwd = 0.7)

####
#### STOCKS
####

iv.cbi3 <- plm(D_logstockvalue ~ LD_logstockvalue + BCsys * cbi + 
                 D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
               + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis   | LD_logstockvalue +
                 D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
               + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis  
               +  L_BCsys * L_CBI_IV_PCA  , 
               data = bc.full, model = "random" )

covDK.st <- vcovSCC(iv.cbi3)
covDK2.st <- sqrt(diag(covDK.st))
stargazer(iv.cbi3,  se=list(covDK2.st), type = "text", star.cutoffs = c(0.1, 0.05, 0.001), star.char = c("*", "**", "***")) 
cbi.st <-  dhplot(model = iv.cbi3, var1 = "BCsys", var2 = "cbi", int = "BCsys:cbi", vcov = covDK.st, ci=0.90,
                  main = "(D) Stocks" , ylab = "(Change) in Stock Capitalization", xlab = "CBI",yint_lwd = 0.7)


CBI.IVs <- ggarrange(cbi.Un, cbi.cr, cbi.st,
                     ncol = 3, nrow = 1)

### Appendix Figure AF3 ###
annotate_figure(CBI.IVs, top = text_grob("Figure AF3: IV Models for CBI, RE with DK Standard Errors",  size = 17))


#### Appendix Table AT2 ####
stargazer(iv.cbi1, iv.cbi2, iv.cbi3,  se=list(covDK2.un, covDK2.cr , covDK2.st), type = "text",
          star.cutoffs = c(0.1, 0.05, 0.001), star.char = c("*", "**", "***"))

##############################################
########### EMP MANDATE MODELS ###############
##############################################

bchighfcbi <- subset(bcs.dat,  FCBI_dum==1)
bc.h.fcbi <- pdata.frame(bchighfcbi, index = c("country", "year"))

####
#### UNEMP 
####

iv.emp1.b <- plm(D_logunemp_ilo ~ LD_logunemp_ilo + BCsys * EmpM + 
                   D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
                 + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis   | LD_logunemp_ilo +
                   D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
                 + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis   
                 +  L_BCsys * L_EmpM_IV_PCA  , 
                 data = bc.h.fcbi, model = "random" )


covDK.un.E <- vcovSCC(iv.emp1.b)
covDK2.un.E <- sqrt(diag(covDK.un.E))
stargazer(iv.emp1.b,  se=list(covDK2.un.E), type = "text", star.cutoffs = c(0.1, 0.05, 0.001), star.char = c("*", "**", "***")) 
Emp.Un <- dhplot(model = iv.emp1.b, var1 = "BCsys", var2 = "EmpM", int = "BCsys:EmpM", vcov = covDK.un.E, ci=0.9,
                 main = "Unemployment" , ylab = "(Change) in Unemployment", xlab = "Employment Mandate Index",yint_lwd = 0.7)

####
#### CREDIT 
####

iv.emp2.b <- plm(D_logdomcredit ~ LD_logdomcredit + BCsys * EmpM + 
                   D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
                 + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis   | LD_logdomcredit +
                   D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
                 + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis   
                 +  L_BCsys * L_EmpM_IV_PCA  , 
                 data = bc.h.fcbi, model = "random" )

covDK.cr.E <- vcovSCC(iv.emp2.b)
covDK2.cr.E <- sqrt(diag(covDK.cr.E))
stargazer(iv.emp2.b,  se=list(covDK2.cr.E), type = "text", star.cutoffs = c(0.1, 0.05, 0.001), star.char = c("*", "**", "***")) 
Emp.Cr <- dhplot(model = iv.emp2.b, var1 = "BCsys", var2 = "EmpM", int = "BCsys:EmpM", vcov = covDK.cr.E, ci=0.9,
                 main = "Credit/GDP" , ylab = "(Change) in Credit/GDP", xlab = "Employment Mandate Index",yint_lwd = 0.7)

####
#### STOCKS
####

iv.emp3.b <- plm(D_logstockvalue ~ LD_logstockvalue + BCsys * EmpM + 
                   D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
                 + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis   | LD_logstockvalue +
                   D_KOFEcGIdf + D_fiscal + L_logcpi_o + L_eldem + L_loggdppc + D_range + L_range + L_EquityBankRatio  
                 + L_KOFEcGIdf + logevol_Reg + D_log_imfgdp2 + L_CrBubble5dum + uscrisis   
                 +  L_BCsys * L_EmpM_IV_PCA  , 
                 data = bc.h.fcbi, model = "random" )

covDK.st.E <- vcovSCC(iv.emp3.b)
covDK2.st.E <- sqrt(diag(covDK.st.E))
stargazer(iv.emp3.b,  se=list(covDK2.st.E), type = "text", star.cutoffs = c(0.1, 0.05, 0.001), star.char = c("*", "**", "***")) 
Emp.st <- dhplot(model = iv.emp3.b, var1 = "BCsys", var2 = "EmpM", int = "BCsys:EmpM", vcov = covDK.st.E, ci=0.90,
                 main = "Stock Mkt", ylab = "(Change) in Stock Capitalization", xlab = "Employment Mandate Index",yint_lwd = 0.7)


Emp.IVs <- ggarrange(Emp.Un, Emp.Cr, Emp.st,
                     ncol = 3, nrow = 1)

#### Figure AF3 Here ####
annotate_figure(Emp.IVs, top = text_grob("Figure AF4: IV Models for Emp Mandates, Re with DK Standard Errors",  size = 14))


#### Appendix Table AT3 ####
stargazer(iv.emp1.b, iv.emp2.b, iv.emp3.b,  se=list(covDK2.un.E, covDK2.cr.E , covDK2.st.E), type = "text",
          star.cutoffs = c(0.1, 0.05, 0.001), star.char = c("*", "**", "***"))


################
## End Script ##
################

